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ANALYSIS OF MULTI SERVER RETRIAL QUEUEING 
SYSTEM: A MARTINGALE APPROACH AND AN 
ALGORITHM OF SOLUTION 

VYACHESLAV M. ABRAMOV 

Abstract. The paper studies a multiserver retrial queueing system 
with m servers. Arrival process is a point process with strictly station- 
ary and ergodic increments. A customer arriving to the system occupies 
one of the free servers. If upon arrival all servers are busy, then the cus- 
tomer goes to the secondary queue, orbit, and after some random time 
retries more and more to occupy a server. A service time of each cus- 
tomer is exponentially distributed random variable with parameter fi\. 
A time between retrials is exponentially distributed with parameter fi2 
for each customer. Using a martingale approach the paper provides an 
analysis of this system. The paper establishes the stability condition and 
studies a behavior of the limiting queue-length distributions as /12 in- 
creases to infinity. As fj,2 — > 00, the paper also proves the convergence of 
appropriate queue-length distributions to those of the associated 'usual' 
multiserver queueing system without retrials. An algorithm for numer- 
ical solution of the equations, associated with the limiting queue-length 
distribution of retrial systems, is provided. 

Keywords: Multiserver retrial queues, Queue-length distribution, Sto- 
chastic calculus, Martingales and semimartingales 

AMS 2000 Subject classifications. 60K25, 60H30. 



1. Introduction, description of the model, review of the literature 

and motivation 

We study a multiserver retrial queueing system having the following struc- 
ture. 

• The arrival process A(t) is a point process, the increments of which form 
a strictly stationary and ergodic sequence of random variables. 

• There are m servers, and an arriving customer occupies one of free 
servers. 

• If upon arrival all servers are busy, then the customer goes to the sec- 
ondary queue, orbit, and after some random time retries more and more to 
occupy a server. 

• A service time of each customer is exponentially distributed random 
variable with parameter n\. 
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• A time between retrials is exponentially distributed with parameter \X2 
for each customer in the orbit. 

Using a martingale approach the paper provides an analysis of this sys- 
tem. The paper establishes the stability condition and studies a behavior 
of the limiting queue-length distributions as ^2 increases to infinity. As 
Hi — > 00, the paper also proves the convergence of appropriate queue- length 
distributions to those of the associated 'usual' multiserver queueing sys- 
tem A/M/m/oo (without retrials), where the first parameter A in the first 
position of the notation denotes the arrival point process A(t). In the follow- 
ing, by 'usual' multiserver queueing system we mean the abovementioned 
A/M/m/oo queueing system. Our asymptotic results can be applied to var- 
ious problems associated with multiserver retrial queues. For example, they 
can be used in performance analysis of real communication systems. 

Analysis of multiserver retrial queueing systems is very hard. For the 
M/M/m retrial queueing systems, analytic results for the stationary prob- 
abilities are not simple even in the case of m = 2. The results associated 
with numerical analysis have been obtained in a large number of papers (see, 
e.g. Anisimov, and Artalejo [4], Artalejo, and Pozo [6], Falin [16], Neuts, 
and Rao [42], Stepanov [48], Wilkinson [51] and others). The methods of 
these papers are based on truncation of the state space for the stationary 
probabilities and constructing auxiliary models helping to approximate the 
initial system (see the review of Artalejo, and Falin [5] as well as the book 
of Falin, and Templeton [17] for details). 

In the present paper we study a non-Markovian retrial queueing system, 
the input process of which is a point process with strictly stationary and er- 
godic increments. By point process with strictly stationary and ergodic 
increments we mean the following. Let {£i}i>i be a strictly stationary 
and ergodic sequence of positive random variables, and let x n = Y17=i 
(xq = 0) be the corresponding sequence of points. Then, the process 
X(t) = Yli^i I{ x i ^ where denotes the indicator of set S, is called a 
point process with strictly stationary and ergodic increments. If {£i}i>i is a 
sequence of independent identically distributed random variables, then X(t) 
is called a point process with independent identically distributed increments 
or renewal process. 

Let E£i = A -1 . Then the assumption that A(t) is a point process with 
strictly stationary and ergodic increments means that 



Along with asymptotic behavior of the limiting queue-length distribution 
as parameter \xi increases to infinity, the paper proves continuity of the 
limiting queue-length distributions in the sense of convergence of functional 
of the queue-length distribution to that of the 'usual' A/M/m/oo queue. 

The significance of our results is motivated as follows. In the queueing lit- 
erature the multiserver M/GI/m/0 and GI/M/m/0 loss queueing systems 
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are the systems of special attention. Both these queueing systems are known 
as a good model for telephone systems and has been an object of investiga- 
tions during many decades. The earliest investigations of these systems were 
due to Palm [43] and later due to Pollatzek [44], Cohen [12], Sevastyanov 
[46], Takacs [49], Iglehart, and Whitt [23]- [25] and others. Recently, the in- 
creasing attention has been to non-stationary multiserver loss systems (e.g. 
Davis, Massey, and Whitt [15], Massey, and Whitt [40]) and multiserver 
systems with multiple customers classes (e.g. Cohen [13], Gail, Hantler, and 
Taylor [18], [19], Righter [46] and others). Both these directions for multi- 
server queueing system are closely related to multiserver queueing systems 
with retrials and abandonments, which recently have been intensively stud- 
ied in the literature in a framework of analysis of call centers (see Garnett, 
Mandelbaum, and Reiman [21], Gans, Koole, and Mandelbaum [20], Grier 
et al. [22], Koole, and Mandelbaum [30], Mandelbaum et al. [36]- [38] and 
others). In a framework of mentioned models of queues, the place of the 
model considered in this paper is clear. Our assumption that input process 
is the process with strictly stationary and ergodic increments is more gen- 
eral than those considered earlier. Even multiserver retrial queueing models 
with recurrent input are very difficult to analysis. The explicit results for the 
stationary distributions of these systems are unknown. The known results 
related to Markovian multiserver retrial queueing models are not sufficient, 
since the real models arising in practice not always can be good approxi- 
mated by Markovian models. Note, that Choi, Chang, and Kim [11] applied 
a not standard MAP\, MAP^jMlm retrial model to cellular networks. 

For such general non-Markovian models as the model considered in the 
paper only the stability results are an object of investigation in the literature 
(see e.g. Altman, and Borovkov [3] and references therein). The most 
relevant model is a model including both retrials and abandonments, as 
it has been considered in the aforementioned papers, associated with call 
centers. We would like to point out that such more extended model can also 
be studied by development of the method of this work. We hope that this 
will be done in the future. 

The analysis of the present paper is based on martingale approach. Nowa- 
days the martingale approach, associated with analysis of different queue- 
ing systems and network, is familiar. Among the well-known general text- 
books on martingale theory such as Jacod, and Shiryayev [26], Karatzas, 
and Shreve [27], Liptser, and Shiryaev [34], [35], Revuz, and Yor [45], there 
are special textbooks on martingale theory associated with point processes 
and queues and networks such as Bremaud [10], Baccelli, and Bremaud [7], 
Whitt [50] and others. Also there is a large number of papers studying dif- 
ferent queueing systems and networks with the aid of stochastic calculus. 
Traditionally, the martingale methods are used to provide weak convergence 
results and diffusion approximations, and the majority of papers establish 
such type of results (e.g. Abramov [1], Kogan, and Liptser [28], Kogan, 
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Liptser, and Shenfild [29], Krichagina [31], Krichagina, Liptser, and Puhal- 
skii [32], Krylov, and Liptser [33], Mandelbaum, and Pats [39], Williams [52] 
and others). 

In some recent papers the martingale theory is used for analysis of point 
processes and queue-length characteristics of queues and networks also under 
the light traffic conditions (e.g. Abramov [1], [2], Kogan, and Liptser [28], 
Miyazawa [41] and references therein). 

In the present paper, we study a behavior of the queue-length process 
under the light traffic condition for the multiserver retrial queueing system, 
by using the known methods of the theory of martingales. The advantage of 
the martingale approach is that, it provides a deepen analysis of the system 
helping to study a more wide its extension, than the traditional methods. 

The paper is structured as follows. There are 11 sections. The main re- 
sults of this paper are given in Section 8. Theorem l8.1l establishes a property 
of the limiting joint probabilities of the queue-length processes in the main 
queue and orbit as parameter \ii increases to infinity. This property is then 
used in Theorem 18.21 stating on the continuity property of the queue-length 
processes, a convergence of the joint queue-length distribution of the multi- 
server retrial queueing system to that of the one-dimensional queue-length 
distribution of the 'usual' A/M/m/oo queueing system. Section 2 discusses 
the basic equations, which are then used throughout the paper. Section 3 
deduces the Doob-Meyer semimartingale decomposition of the basic equa- 
tions. Section 4 studies normalized queue-length processes and establishes 
the condition for stability (in the sense defined precisely in Theorem 14. lj) . 
Section 5 derives equation for the queue-length distribution given by Theo- 
rem 15.11 Section 6 is devoted to the proof of Theorem 15.11 There are two 
corollaries of Theorem 15. II given in Section 7. Sections 9 and 10 discuss al- 
gorithm for numerical solution of the main system of equations. Specifically, 
Section 10 provides numerical results under conditions of high retrial rate, 
which enable us to discuss the results on convergence obtained in Section 8. 
Concluding remarks are given in Section 11. 



All point processes considered in this paper are assumed to be right- 
continuous having the left-side limits. 

Consider the queue-length process of our retrial system. The number of 
servers, occupied in time t, are denoted Qi(t), and respectively, is the 

number of customers in orbit in time t. The both queue-length processes 
Qi(t) and Q2(t) are assumed to be continuous in 0, Qi(0) = (^2(0) = 0, 
as well as right-continuous having the left-side limits. The following two 
equations describe a dynamic of the queue-length processes Qi(t) and Q2{t). 
The first equation is 



2. Discussion of the basic equations 



(2-1) 
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where tt\ , i = 1, 2, ...,m, are independent Poisson processes with rate 
The second equation is 



where 7r) , % = 1, 2, are independent Poisson processes with rate [i<i- 

Equations (|2~T|) and lf2~2^l can be explained as follows. The term ^2i=i 
I{Qi(s— ) > i} of the integrand of (|2.1|) means the number of occupied 
servers immediately before time s in the main queue. We use the term 
'immediately before' keeping in mind that the point s can be a point of 
possible jump. Then the right-hand side of equation (|2.1jl for Q\(t) + Q2(t) 
includes the number of arrivals until time t minus the number of departures 
given by the term 



The term X/Si ^-{Q"2{ s ~ ) > of the second integrand of (|2.2[) means the 
number of customers in orbit immediately before time s. Obviously, if there 
is no customer in orbit, then the second integrand of ()2.2j) becomes equal to 
0. Next, the term \{Q\{s— ) ^ m} of the second integrand of (|2.2j) means 
that if immediately before time s there is at least one free server in the 
main system, then one of the customers of the orbit queue can occupy the 
server in time s, otherwise the integrand becomes equal to 0. The first 
integral of the right-hand side of ()2.2|) means the number of arrivals to the 
orbit system during the time interval [0,t]. Then the right-hand side of 
equation (|2,2j) for Q2{t) includes the number of arrivals until time t to the 
orbit minus the number of departures from the orbit to the main queue, 
where the mentioned number of arrivals to the orbit is given by the first 
integral, and the mentioned number of departures from the orbit is given by 
the second one. 

Notice, that the equations similar to (|2.1j) and (|2.2|h associated with time- 
dependent Markovian model with abandonments and retrials, has already 
been considered in the literature (e.g. Mandelbaum et al. [36]- [38]). How- 
ever, the main emphasis of these papers was done to the analysis of fluid 
limits and diffusion approximations. 

3. Semimartingale decomposition of the queue-length process 

In this section we provide another representation for the queue-length 
processes by using the Doob-Meyer semimartingale decomposition (see e.g. 
Liptser, and Shiryayev [34], Jacod, and Shiryayev [26]). The compensators 
of the semimartingales associated with point processes will be provided by 
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(2.3) 
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'hat'. For example, the point process A(t) is a semimartingale, and A(t) is 
its compensator. The Doob-Meyer semimartingale decomposition for some 
semimartingale X will be written as X = X+Mx, where Mx is the notation 
for the local square-integrable martingale. For example, the semimartingale 
decomposition of A(t) is written as A(t) = A(t) + M A {t). Along with the 
notation Mx sometimes it is also used Mj(t) or Mjj(i). In such cases the 
sense of the local square integrable martingales Mj(t) or M{j(t) is specially 
explained. 

It is assumed in the paper that all point processes are adapted with respect 
to the filtration given on stochastic basis {O, JP, F = (J^i)t>o, P}- 
Let us start from equation (|2.1|) . As semimartingales, the processes A(t) 

and C(t) = Jq Y^JILi I{Qi( s— ) — *'}d7r| ^ (s) are represented as 

(3.1) A(t) = A(t)+M A (t), 

(3.2) C(t) = C(t)+M c (t), 
The compensator C(t) has the representation 

(3.3) C(t) = in [ Qi(s)ds 

Jo 

(for details see Dellacherie [14], Liptser, and Shiryayev [34], [35], Theorem 
1.6.1). 

By virtue of (|3.1|) . (|3.2|) . and (|3.3j) . for equation i|2.1j) . we have 

(3.4) Q x {t) + Q 2 (t) = A(t) + M A (t) - m f Q l {s)ds - M c {t). 

Jo 

Denoting Mi{t) = M A {t) - M c {t) we obtain 

(3.5) Qx{t) + Q 2 (t) = A(t) - m I Qi(s)ds + M x {t). 

Jo 

Let us now consider equation (|2.2[) . For the associated arrival process 
D\{t) = J I{Qi(s— ) = m}dA(s) we have the following: 

(3.6) D 1 (t) = D l (t) + M Dl {t), 
where 

(3.7) 5i(t)= / I{Qi(«-) = m}<L4(s) 

Jo 

and 

(3.8) M I3l (i)= /"*I{Q 1 ( S -)=m}dM A (s). 

JO 
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For the associated departure process D2(t) = J I{Qi(s— ) ^ m} Yli^i 

(2) 

1{Q2\S— ) > i}dir^ (s) we have the following: 

(3.9) D 2 (t) =D 2 (t)+M D2 (t), 
where 

(3.10) D 2 (t) = » 2 I I{Qi(s) ± m}Q 2 (s)ds 

Jo 

(see Dellacherie [14], Liptser, and Shiryayev [34], [35], Theorem 1.6.1). 
Then (12.21) can be rewritten as follows: 



(3.11) 



Q 2 (t)= [ I{Q 1 (s-) = m}dA(s) 
Jo 

-fi 2 f I{Q l (s)^m}Q 2 (s)ds + M 2 (t), 
Jo 



where 

(3.12) M 2 (t) = M Dl {t) - M D2 {t). 

4. Normalized queue-length processes and condition for the 

stability 

In this section we study the normalized queue-length processes 

(4.1) qk(t ) = 9M : fc = 1,2; t > 0, 

and its asymptotic properties as t — > oo. In the following the small letters 
will stand for normalized processes. The notation for normalized processes 
corresponds to the notation of original processes given by capital letters. 
For example, a(t) = t^ 1 A(t); mn 2 (t) = t~ l MD 2 (t) and so on. 

Following this comment, the equations associated with the queue-length 
processes (|3.5j) and (J3.ll)) can be written 

(4.2) q x (t) +q 2 (t) =a(t) - ^ f s qi (s)ds + mi(t), 

1 Jo 

and 



(4.3) 



q 2 {t) = - f I{Q 1 (s-) = m}d[sa(s)] 
T Jo 

I{Qi(s) / m}sq 2 (s)ds + m 2 (t), 



t 



Let us now study these two equations (|4.2|) and (|4.3I) as t — > oo. More 
accurately, let us find Plim of qi(t) and q 2 (t) as t — > oo, where Plim denotes 
the limit in probability. 

Show that 

(4.4) Plimmi(t)=0. 

t^oo 
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Indeed, because of m\(t) = m A (t) — mc(t), we have 



lim mi(t) < P lim m^t) + P lim mc(t) . 



t— >oo 



(4.5) 

Applying the Lenglart-Rebolledo inequality we obtain: 

sm A {s) 

0<s<t 



F{\m A (t)\ > 5} < P{ sup 



t 



> 6 



} 



(4.6) 



{ sup \M A (s)\ > St\ 

0<s<t > 



<^+P{A(t)>et 2 } 



P- 



+ 



t 



> et 



The both terms of the right-hand side vanish, as e is taken sufficiently small 



and t increases to infinity such that et — » oo. That is Plimt_ +0O 
Taking into account that 

ft m m 

(4.7) 



\m A (t)\ =0. 



nt llli III' 

/ £l{gi(«-) > *}d^ {1) (s) < $>«(*), 
i=i i=l 



by virtue of the Lenglart-Rebolledo inequality we have: 

smc(s) 



P{\mc(t)\ >5} < P{ 



sup 

0<s<t 



(4.8) 



= P{ sup U/ C (s)| > <5tj 

m 

<^{^>^ >et2 } 

i=l 

+p{J^« (1) w>4 

| m c(*)| = 0- Thus, 



e 

P 



i=l 



As earlier (see reference (jUSJ)), now we obtain Plim^ — > rC<) 
it is shown that Plim^oo m\ (t) = 0. 

Analogously to the above, m,2(t) = m_o 1 (t) — mu 2 (t). Therefore, similarly 
to (HSl) 



(4.9) 



lim m2(t) <P lim mc^t) +P lim mo 2 (i) . 



t— >oo 



t^oo 



t—fOO 



Notice (see (|3.8|l ). that \mD 1 (t)\ < \m A (t)\ for all t > 0. Therefore, 
(4.10) P lim |m Dl (i)| < P hm \m A (t)\ = 0. 

However, both Plim^oo |m/) 2 (i)| = and Plim^oo 1 777,2 (*) I = can only 
be true under the additional condition Plim^ 00 a(t) < [X\m. Recall that 
according to (jl.ll) Plim^ 00 a(t) = A. Therefore, the above condition is 
A < fiim. 
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In order to prove Plim^oo |m2(i)| = and Plim^oo |m£) 2 (t)| = under 
the abovementioned additional condition A < [i\m let us now study equation 

dOJ). 

Notice first that from the fact that Plimt_ i00 \mji(t)\ = we have 

(4.11) P lim a(t) = P lim a(t) = A, 

since 

(4.12) P lim a(t) = P lim a(t) - P lim m A {t). 

t^oo t^oo t^oo 

Next, if limt_»oo i -1 Jo P{Qi (s) = m}ds = 1, then by the Lebesgue domi- 
nated convergence 

(4.13) E lim - f I{Qi(a) = m}d[sa(s)} = A, 

t^QO t J 

and 

(4.14) lim Eq 2 (t) = lim Eo(t) = A. 

(gllg) means that if lirn^oo / ' P{<5i(s) = "i} ds = 1, then EQ 2 (t) in- 
creases to infinity, as t — > oo. 

Let us now assume, that lim^oo^ 1 / * P{Qi (s) = m}ds = p < 1. Then, 
it is not difficult to prove that only Plim^oo q 2 (t) = must satisfy (|4.2[l . 

Indeed, assume that P linn,—^ q 2 (t) > 0. Then, taking into account that 
lim^oo t^ 1 J Q ¥{Qi(s) 7^ m}ds = 1 — p > 0, for large i we obtain 

(4.15) I l{Qx{s) ^ m}sq 2 {s)ds = 0{t 2 ). 
Jo 

This means that 



oo, 



(4.16) P lim ^ I I{Qi(s) ^m}sq 2 (s)ds 

and therefore lim^oo Eq 2 (t) = — oo. This contradicts to the fact that 
Eq 2 (t) > for all t > 0, and therefore, only Plim^oo <?2( s ) = is a pos- 
sible limit in probability, satisfying (|4.3|) for some unique value p = p* = 
limt-,00 t- 1 J * P{Qi(s) = m}ds. 

Thus, we proved that if lirm ; _ i . 00 J„ P{Qi(s) = m}ds = p < 1, then 
Plim^oo 92(5) = 0. In fact, taking into account that Q 2 (t) is a cadlag 
process having with probability 1 a finite number of jumps in any finite 
interval, from the above contradiction we have 

1 r* 

(4.17) lim - / EQ 2 (s)ds < 00, 

t^OQ t J Q 

and therefore, 

(4.18) P lim - f I{Q 2 (s) < oo}ds = lim - f F{Q 2 (s) < oo}ds = 1. 

t^oo t J t^co t Jo 
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Now, using the Lenglart-Rebolledo inequality for large t we obtain 



P{|m D2 (i)|><5}<p{ 



sup 

o<s<t 



sm D2 {s) 



t 



> 6 



} 



"{ sup \M D2 (s)\ > St} 

0<S<t 



(4.19) 



^^2 +¥ {f o HQlis) + m}sq{s)ds > et 2 } 
< ^ + P{^ sq{s)ds > et 2 } 

=i +¥ {ii! sq{s)ds>et y 



Therefore, under the assumption that linit_ >00 t 1 Jjj P{Qi(s) = m}ds = p < 
1, we obtain 

(4.20) 



P lim \m D2 {i) \ = 0. 



t— >oo 



Hence, together with (|4.10j) . under the assumption that lim t ^oo t 1 J c 
P{Qi(s) = m}ds = p < 1 we have 

(4.21) P lim \m 2 (t)\ = 0. 

Let us return to relation (|4.18f) under the condition A < [i\m. We have 



(4.22) 



lim - / VP{Q 2 (s)=i}ds = l, 

:^0O t In * — ' 

J u 1=0 



and, since Qi(s) can take values 0, 1, . . . , m only, we have 

l-t m 

(4.23) 



lim - / VPWi( S )=I}d S = l. 



Therefore, under the condition A < Him, 

m oo ct 



VVlim- P{Q!( S ) = Z,Q 2 ( S ) = fc}d s 

Z ' Z ' i^OO t In. 

1=0 k=0 J{) 

(4 24) 1 ft m °° 

= lim- / VVP{Q 1 ( S ) = /,Q 2 (s) = A : }d S 

J U /— n z — n 



l. 



Thus, we have the following theorem. 

Theorem 4.1. Under the condition A < mm there exist 

1 /■* 

lim - / P{Qi(s) = i,Q 2 (s) = fc}ds, 
I = 0, 1, . . . , m; k = 0, 1, . . . , 
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satisfying (|4.24[) . 



5. Analysis of the limiting queue-length distributions 

In the rest of the paper it is assumed that condition A < is fulfilled, 
and therefore the system is stable in the sense of Theorem 14.11 Notice, 
that the statement of Theorem 14.11 does not mean existence of the limiting 
stationary probabilities as t — > oo. For example, if increments of the point 
process A(t) are lattice, then the stationary probabilities do not exist. 

In this case one can speak only about appropriate fractions of time in two- 
dimensional states (i,j) associated with the queue-length processes Qi(t) 
and Q2(t). In general, we speak about 



rather than lim t _ too ¥{Qi(t) = i, Q2W = j}- However, if the increments 
of the point process A(t) are independent, identically distributed and non- 
lattice, then there exist the limiting stationary probabilities lim^oo ¥{Qi(t) = 
i, Qzit) = j} coinciding with (|5.1|) . Indeed, in this case the process Q(t) = 
Qi(t) + Q2{t) has a structure of regeneration process, and then the proof 
of this fact follows by a slight extension of arguments given in the proofs of 
Theorem 5 on p. 173 and Theorem 22 on p. 157 of Borovkov [9]. 

Let us introduce the processes 



(5.2) I iJ (t) = I{Q 1 (t) = i n Q 2 (t)=j}, i = 0,l,...,m; j = 0,l,... , 



assuming that I_ij(t) = and Ij ) _i(t) = 0. 

The jump of a point process is denoted by adding A. For example, AA(t) 

is a jump of A(t), A7r[^(t) is a jump of the fcth Poisson process with rate 
Hi etc. 

Let us denote: 



(5.1) 




(5.3) 




k=l 



and 



3 



(5.4) 




k=l 
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Taking into account that the jumps of all the processes A(t), ir^ (t) and 

7iy (t) are disjoint (k = 1,2, ...,m; Z = 1,2, ...), we have the following equa- 
tions: 

HQi(t-) + AQ 1 (t) = i n Q 2 (t) = Q 2 (t-) = j} 



(5.5) 



(5.6) 



Ii- ld (t-)AA(t) + (t-)An^(t) 



+ j w (t-)[i - AA(t)][i - An^(t)][i - Anf (t)], 

i = 0, 1, m — 1; 

I{Qi(«-) + AQi(t) = i n Q 2 (t-) ^ Q 2 (t) = j} 
= j i _i J -+i(*-)Ang 1 (*), 

z = 0, 1, m — 1; 



I{Qi(*-) + AQi(*) = m n Q 2 (t) = Q 2 (t-) = j} 
(5.7) =I m _i, 3 (t-)Ai(t) 

+ / mJ H[i-A#)][i-AnW(t)] ; 



(5.8) 



I{Qi(*-) + AQi(t) = m n Q 2 (t-) ^ Q 2 (t) = j} 
= I mJ -i(t-)AA(t). 



Then, 



(5.9) 



AI itj (t) = I{Qi(*-) + AQi(i) = i n Q 2 (t) = Q 2 (t-) = j} 
+ I{Qi(t-) + AQi(t)=t n Q 2 (t-) ^ Q 2 (t) = j} 

-IiAt-), 

i = 0,1, ...,m; j > 0. 



Since 



(5.10) 



j;A/y( S )=/y(t)-Jy(0), 



we have the following. 



(5.11) 
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For i = 0, 1, m — 1, 

hj{t) = I itj (0) + f [Ii-ij(s-) - / M -( S -)]cL4( S ) 
Jo 

- jf* j i)i ( s -)dnf ) ( s ) - jT J iJ ( s -)dnf ( a ) 
+ / / i+1J ( s -)dn|J 1 ( s ) 

JO 

+ / J i -i J+ i( s -)dnfi( s ). 

" 

In turn, for i = m we have 

ViW = + / [/m-i,j(s-) - I mj -(a-)](LA(s) 



(5.12) - /' < / m)j ( s _) d n«( s ) + /'*J mJ _ 1 ( s -)dA( S ) 

j •/ 

+ / / m - 1J+ i( s -)dnfi( s ). 

•/ 

Using the Doob-Meyer semimartingale decomposition, from (|5.11j) and 
(|5.12j) we obtain the following equations. For i = 0, 1, m — 1 

I id (t) = 1^(0) + / [li-ij(a-) - Ii,j(s-)}dA(s) 



-fill / Ii,j(s)ds - fj, 2 j Ii,j(s)ds 
(5.13) • /o t 1/0 

+ + 1) / / i+ ij(s)ds 

JO 

+ /i 2 (j + l) / Ji-ij+iWds + Mijft), 
Jo 

where the local square integrable martingale Mjj(t) has the representation 

MiA*) = f [Ji-ij(a-) - iij(a-)]d[A(«) - l(s)] 
Jo 

i M -( s -)d[n l (1) ( s )- / i 1 i S ] 



3 

t 



(5.14) - jf / i)j ( s -)d[nf( s )_ /X2is ] 

+ / / i+y ( s -)d[n« ( s )-(i + i) ms ] 

Jo 

+ / / i -i J -+i( s -)d[ng 1 ( s )-(i + i) Al2S ]. 

« 
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In turn, for i = m we have 

Imj(t) = I m ,j(°) + [ [Im-l,j(s-) ~ I m j{8-)]dA(8) 

Jo 

(5.15) -mm \ I m j(s)ds+ / 7 mJ _i(s-)dA(s) 

J ■/ 

+ /x 2 (j + l) I I m -i,j+i{s)ds + M md (t), 
Jo 

where the local square integrable martingale M m j(t) has the representation 



(5.16) 



i(t)= [ [I m -lj(*-) ~ Imj(s-)]d[A(8) - A(8)] 

Jo 

- I m ,j{s-)d[n$(s) -turns] 
Jo 

+ / I m ,j-i(s-)d[A(s) - A(s)] 
Jo 

+ I I m -i, j+ i{s-)d[Yl%{s)-tJi 2 {j + l)s]. 

J 



Now, we are ready to formulate and prove the theorem. 



Theorem 5.1. For the queue-length processes Qi(t) and Q 2 (t) we have the 
following equations. 



(i) In the case i = 



1 /■* 

lim -E / I{Q 1 (s-) = 0,Q 2 {s-)=j}dA{s) 

t^oo t Jq 

1 i 4 

(5.17) = hi lim - / P{Qi(s) = 1, Q 2 (s) = j}ds 

t^oo t Jq 

i r* 

-ti 2 j}im - / P{Qi(s) = 0,Q 2 (s) = 
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(ii) In the case i = 1,2, m — 1 (m > 2) 

lim fmQ^s-) = i, Q 2 {s-) = j} 

t^oo t Jq 

-I{Q 1 (s-)=i-l,Q 2 (s-)=j}]dA(s) 
1 /"* 

= m(i + 1) lim - / P{Qi(s) = i + 1, Q 2 (s) = j}ds 

t— »oc t Jq 
(5.18) . 1 rt 



fiii lim - / P{Qi(s) =i,Q 2 (s) = j}ds 

t^OO t Jq 

-^j lim- / F{Q 1 (s) = i,Q 2 (s)=j}ds 

t~>OC t Jq 

1 /"* 

+ /x 2 (j + 1) lim - / P{Qi(s) = * - 1, Q 2 (s) = j + l}ds. 

t^OO t Jq 

(Hi) In the case i = m 

lim / [I{Qi(s-) = m,Q 2 (0 =j} 

i-^OO t Jq 

-I{Q 1 (s-)=m-l,Q 2 (s-)=j} 

(5 19) ~ I{Ql(s_) = m ' g2(s_) = j ~ 1}](M(s) 

1 f l 

= fi 2 (j + 1) lim - / P{Qi(s) = m - 1, Q 2 (s) = j + l}ds 

t^OO t Jq 
1 t* 

- mm lim - / P{Qi(s) = m, Q 2 (s) = j}. 

t^OO t Jq 

Here in H5.17j) - ()5.19[) it is put I{Q\(t) = i, Q 2 (t) = j} = if at least one of 
the values i, j is equal to -1. 

6. Proof of Theorem 15.11 

Let us first study the case i = 0. From (|5.13|) we have 

Phm I(/ .(t)_/ .(0)) = -Plim - f ' I j(s-)dA(s) 

+ P lim -/xi / h,j(s)ds 

(Q I) * JO 

1 /■* 

-P lim -/i 2 j / Joj(s)ds 
+ Phm 

The left-hand side of (|6.1f) is equal to zero. Therefore, rewriting (|6.1|) in the 
form = —K\ + K 2 — K3 + K4, let us compute the terms of the right-hand 
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side. Using the Lebesgue theorem on dominated convergence we have 



1 



(6.2) 



Kt = P lim - / I 0J (s-)dA(s) 



lim -E / I j(s-)dA(s). 



Taking into account that A{t)/t and A(t)/t have the same limit in probability 
(see (|4.11(l ). relation (|6.2[) can be finally rewritten as follows: 



(6.3) 



Next, 



(6.4) 



i r* 

K x = lim -E / = 0, Q 2 {s-) = j}dA(s) 

= lim ±E / I{Qi(s-) = 0, Q 2 (s-) = j}d[A(s) - M A (s)} 
= lim ^E f I{Qi(s-) = 0, Q a ( S -) = j}dA(s). 



i r* 

K 2 =¥ lim -m / hj(s)ds 

t^QO t Jq 

1 /** 

= Ml lim - / P{Qi(s = 1, Q 2 (s) = j}ds. 



Similarly, 



(6.5) 



i r* 

P - hm ~(j, 2 j / h,j(s)ds 

t^oo t J 

i r* 

» 2 j lim - / P{Qi(s) = 0,Q 2 (s) = j}ds. 



Notice, that if j = then iT 2 = 0. Next, 

Mi,j(t) 



K 4 = P lim 

t— >oo 



(6.6) 



t^oo 



< P lim (\m A (t) \ + 



^\t)-^t 



+ 



+ 



(2) 



o. 
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Thus, for i = we have the following 

1 I* 

lim -E / I{Qi s- =0,Q 2 (s- =j}dA{s) 

t^oo t Jq 

(6.7) = Ml lim 1 / P{Qi(s) = l,Q 2 (s) =j}ds 

1 /** 

- M2 j lim - / P{Qi(a) = 0, Q 2 (s) = j}ds. 

()5.17j) follows. 

Let us consider now the case 1 < i < to — 1 (to > 2). We have the 
following equation: 

P lim l(li,j(t) - I id (pj) = P lim I [ Ii_ ld (s-)dA(s) 
t— >oo (. \ / t^oo t Jq 

1 f 4 

-P lim - / Iij(s-)dA(s) 

t^CO t Jq 

1 /■< 

+ P lim 7 /ii(t + l) / / i+ i J (s)ds 
* Jo 

1 /"* 

(6-8) -P lim -fj, 2 j / /t,-(a)ds 

t^OO t Jq 
1 /** 

— P lim -fj,\i / 7j,-(s)ds 

t^oo t Jq 

1 /■* 

+ P lim -M2(j + 1) / Ii-ij+i(a)ds 

i^OO t Jq 
M i i(t) 

+ P lim l,n ' , 

t^oo t 

The term of the left-hand side of (|6.8f) is equal to zero. Therefore rewriting 
(|B~H)) as = Ki - K 2 + K 3 - K± - K 5 + K 6 + K 7 we have the following. 
Similarly to ()6.3|) 

(6.9) K x = lim -E /* I{Qi(*-) = i - 1, Q 2 (s-) = j}dA(s), 

t^oo t Jq 

and 

1 f* 

(6.10) K 2 = lim -E / I{Qi(s-) = i, Q 2 (s-) = j}cL4(s). 

t^oo t Jq 

Similarly to (fO|) and (fB3f) 

1 /■* 

(6.11) K 3 = m(i + l) lim - / P{Qi(s) = i + l,Q 2 (s)= j}ds, 

i r* 

(6.12) iT 4 = M 2 i lim - / F{Q 1 (s) = i,Q 2 (s)=j}ds, 

t^oo t Jq 
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1 /"* 

(6.13) K b = mi lim - / P{Qi(a) = t, Q 2 (s) = j}da, 

1 /■* 

(6.14) K 6 = M2 (j + 1) lim - / P{Qi(s) = i - 1, Q 2 (s) = j + l}ds, 

t^oo t Jq 

and similarly to (|6.6|) 

(6.15) K 7 = 0. 
Thus, for 1 < i < m — 1 (m > 2) we have 

lim [ \t{Qi(8-) = i,Q 2 (8-)=j} 

t^oo t Jq 

- l{Q x {a-) = i - 1, Q 2 (s-) = i}]d^(s) 

1 /"* 

= + 1) lim - / P{Qi(s) = » + 1, Q 2 (s) = j}ds 

(6.16) i ft 

- mi lim - / P{Qi(s) = i, Q 2 (a) = j}ds 

1 /"* 

- M2J lim - / P{Qi(s) = i, g 2 («) = J'jds 

1 /"* 

+ /x 2 (j + 1) lim - / P{Qi s = % - 1, Q 2 (s = j + l}ds. 

(5.18) follows. 

Now, consider the last case i = m. We have 

P lim t(W*)-W0)) 
1 /•* 

= P lim - / / m _ lij S - di s) 
-P lim i / I m>j (s-)dA(s) 

t^OO t Jq 
1 /"* 

(6.17) + Plim-/ / mj _i(s-)dA(s) 

t Jo 

1 /•* 

- /iimP hm - / [/ m j(s) - I m ,j-i(s)]ds 
t-*oo t Jq 

1 /"* 

+ /x 2 0' + 1)P lim - / / m _i ji+ i(s)ds 
+ PUm 

t^oo t 

The term of the left-hand side of ()6.17j) is equal to zero. Therefore rewriting 
(|6.17j) as = K\ — K 2 + K% — K4 + K 5 + Kq analogously to the above cases 
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we have the following: 

(6.18) K x = lim -E / I{Q 1 (s-)=m-l,Q 2 (s-)=j}dA(s), 

t^oo t Jq 

i r* 

(6.19) K 2 = lim -E / I{Qi(s-) = m, Q 2 (s-) = j}dA(s), 

t^OO t Jq 

i r* 

(6.20) K 3 = lim -E / 1{Q X s- = m, Q 2 (s- = j - l}dA(s), 

t^oo t Jq 

i r* 

(6.21) K 4 = mm lim - / P{Qi(s) = m, Q 2 (s) = j}, 

1 /■* 

(6.22) K 5 = M2 (j + 1) lim - / P{Qi(s) = m - 1, Q 2 (s) = j + l}ds, 

(6.23) K 6 = 0. 
Thus, for i = m, we have the following: 

i r* 

lim -E / [I{Qi(s-)=m,Q 2 (s-)=j} 

<->oo t Jq 

(6.24) 



-l{Q 1 (s-) = m-l,Q 2 (s-)=j} 

- I{Qi(s-) = m, Q 2 (s-) = j - l}]dA(a) 



= M2(j + 1) lim 7 / P{Qi(s) = m - 1, Q 2 (s) = j + l}ds 

1 /■* 

- [i x m lim - / P{Qi(s) = m, Q 2 (s) = j}. 

t^OO t Jq 

(|5.19|) follows. Theorem 15. II is proved. 

7. Special cases 

Two corollaries of Theorem 15.11 are provided below. The first corollary 
is related to the special case when the process A(t) is Poisson. This case 
is well-known and can be found in Chapter 2 of the book of Falin, and 
Templeton [17]. The second corollary is related to the case of the 'usual' 
A/M/m/oo queue. 

Corollary 7.1. If A(t) is a Poisson processes with rate X, then we have the 
following system of equations, 
(i) In the case i = 

(7.1) XPqj = HiPij - ^jPoj- 
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(ii) In the case i = 1,2, m — 1 (m > 2) 

= fil(i + - Ml^-Pij - M2i-Pjj + M2(j + 

(Hi) In the case i = m 

^ ^{Pm,j — Ptn-lJ ~ PmJ-l) 

= + l)-Fm-l,j+l - HifnPm.j 

Here in l|7.1|) - ([7.3|) we use the notation Pi j = lim^oo P{Qi(i) = i, Q%{t) = 
3}. 

Proof. The proof of Corollary 17.11 follows easily from the statement of The- 
orem 15.11 Indeed, taking into account Ijfi.MJI and the fact, that when the 
process A{t) is Poisson with rate A, then we have A(t) = Xt. Finally, the 
result follows by taking into account the existence of the limiting stationary 
in time probabilities, that is 



Pij = lim F{Q 1 (t)=i,Q 2 (t)=j} 

t — ►oo 

( 7 - 4 ) = lim \ f P{QiO) = i, Q 2 (s) = j}ds, 

i = 0,1 5 ... ,m; j = 0,1, .... 



□ 



Corollary 7.2. LetQ(t) denote the queue-length process for the multiserver 
queueing system A/M/l/oo. Then we have the following system of equa- 
tions. 

(i) In the case i = 



(7.5) 



1 /"* ~ 

lim -E / I{Q(s-) = 0}dA(s) 
-»oo i ,/ 

/xi lim - / ¥{Q(s) = l}ds. 
t^oo t J 



(ii) In the case i = 1,2, m — 1 (m>2) 
1 /■* ~ 

lim -E / = z} - I{Q(s-) = « - 1} dA(s) 

1 r* ~ 

(7.6) = ^(i + 1) lim - / P{Q(s) = ? + l}cfe 

1 /"* ~ 
- mi lim - / ¥{Q(s) = i}ds. 

t~>co t J 
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(Hi) In the case i > m 



lim -E [ [I{Q(s-) =i}- I{Q(s-) = i - l}]dA(s) 

^oo t J 

1 /■* 

/iim lim - / [P{Q(s) = j + 1} - P{Q(s) = i}) ds. 

t^OO t J Q 



(7.7) "~' 



Proof. In order to prove Corollary 17.21 notice that the queue-length process 
satisfies the following equation 

ft °° 

(7.8) Q(t) = 



r ' — ' 

/ ^I{Q( S -)>,}dvr 4 (1) ( S ), 
J ° i=l 



where, as earlier, {7r^ } is a sequence of independent Poisson processes with 
rate [i\ . Then the proof of Corollary 17.21 is analogous to that of Theorem 
15.11 and based on the following equations: 

Ali(t) = Ii-i(t-)AA(t) + I i+1 (t-)AD£\(i) 

+ A/ f (t-)[l - AA(i)][l - AnfV)] - /;(*-), 

where = = i}. □ 



8. Asymptotic analysis of the system as ^2 increases to infinity 

In this section we study a behavior of the system as fj,% increases to infinity. 
Specifically we solve the following problem. 
• How behave 

lim - fnQi(s) = i,Q 2 {s)=j}&s 

when i < m, j > 1? 

The answer to the above question is given by the following theorem. 
Theorem 8.1. As [12 — ► oo ; then for all i = 0,1, m — 1 



00 -\ ft 

- 1 ) / V{Qi{s)=i,Q 2 {s)=j}ds S ) =0(4- 

j=l J o 

Proof. Notice first, that by virtue of (|4.17j) 

m 00 1 /"* 

( 8 - 2 ) EE-K/i™ 7 / F WiW=^2(s) = i}d S )<a). 

i=0i=l 1/0 
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Let us now start from (|5.17j) . Dividing this equation to large parameter /x 2 
(j > 1), we obtain 



i r* 

lim - / P{Qi(s) =0,Q 2 (s) = j}ds 
t Jo 

(8.3) = ^ lim i / P{Qi(a) = l,Q 2 (s) = j}ds 

-J-lim ±E / I{Q 1 ( S -) = 0,Q 2 (s-)=i}dA(s) 
Therefore, from 1)8.3(1 we obtain 

i r* 

j lim - / P{Qi(s) = 0,Q 2 (s) = j}ds 

(8.4) 1/0 ( 

< ^ lim 1 / P{Q 1 ( a ) = 1,Q 2 ( S ) =j}ds > 

with an absolute constant Cq satisfying Co < f-ii- Therefore, 



i.5) 



OO 1 ft 

^ (j lim - / P{Qi(s) = 0, Q 2 ( S ) = j}ds^ 
j=i ^° 

^— E(^ lim 7 / P{QiW = i,Q 2 ( s )=i}d s ). 



.7 = 

Notice now that because of as /z 2 ~ *• 00 > the expressions 

i r* 

(8.6) lim - / P{Qi(s) = i, Q 2 (s) = j}ds, 

t Jo 

and 

(8.7) lim 1e / I{Qi(s-) =i,Q 2 (s-)=i}<L4(s) 

are of the same order. Then considering equation ()5.18|) divided as earlier to 
large parameter // 2 (j > 1), with the aid of induction by the same manner 
we obtain 

00 / 1 /•* 

h(s) = i,Q 2 (s) = j}ds' 



OO 1 /"i 

^(7 lim - / P{Q 

^ V t->oo t Jo 

(8-8) 

<-V Jlim - / F{Q 1 (s) = i + l,Q 2 (s)=j}ds 

where Cj is an absolute constant, i = 1,2, ...,m— 1 (m > 2). The statement 
follows. □ 

Theorem 8.1 helps us to establish the continuity theorem. Denoting 

1 /"* 

(8.9) Ji = lim - / P{Qi(s) = i, Q 2 (s) = J'}ds, 
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and 

1 /■* ~ 

(8.10) J 2 = lim - / P{Q(s) = i + j}ds, 

t^OO t Jq 

we have the following. 

Theorem 8.2. Let \i 2 — > oo. T/ien in i/te cases (7) j = and (2) i = m 

the difference between J\ and J 2 is o(l). 

Proof. It is clear that the difference between J\ and J 2 should be studied 
only in the two cases mentioned in the theorem, since in other cases, as 
Hi — > 00, J\ tends to while J 2 remains positive in general. 

In the case j = we have the following. When i = 0, ()5.17|) coincides 
with (|7.5|) . When i = 1, 2, . . . , m — 1, m > 2, from (|5.18|) we have 

1 /"* 

lim -E / [I{Qi(s-) = i,Q 2 («-) = 0} 

- I{Qi(s-) = i - 1, Q 2 (s-) = 0}]d,4(s) 

1 /"* 

(8 n) = Mi(* + 1) ^ 1 j = i + 1, = 0}ds 

1 /"* 

- pit lim - / P{Qi(s) = i, Q 2 (s) = 0}ds 

1 /■* 

+ a 2 lim - / P{Qi(s) = i - 1, Q 2 (s) = l}ds. 
* Jo 

According to Theorem 18.11 the last term of the right hand-side vanishes as 
(JL2 — > 00. Therefore the limiting relation, not containing this last term, 
coincides with corresponding relation of (|7.6jl . 

In turn, in the case i = m and j > 1 from (|5.19j) we have 



1 /■* 

lim -E / [I{Q 1 (s-)=m,Q 2 (s- 



j} 



- I{Qi(s-) = m, Q 2 («-) = j - l}]dA(a 
U2) 1 f* 



/i 2 (j + 1) lim J / P{Qi(s) = m - 1, Q 2 (s) = j + l}ds 

t^oo t Jq 

-mm lim - f F{Qi(s) =m,Q 2 (s)=j}ds + o(—). 

t^oo t Jq V/i 2 / 

Denoting Q(t) = m + Q 2 (i), then (ETHl can be rewritten 
1 /"* 

lim -E / [I{Q(s-) = m + j} - I{Q(s-) = m + j - 1} ]dA(i 

t^QO t Jq 

1 /"* 

(8.13) = M2 (j + 1) lim - / P{Qi(s) = m - 1, Q 2 (s) = j + l}ds 

t^OO t Jq 

- aim lim - f P{Q(s) = m + j}ds + o(-) . 

t^CotjQ V/X 2 / 
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Now, comparison with (|7.5jl - (|7.7|) and the normalization condition enables 
us to conclude that 

ft 



i r 

lim fi 2 (j + 1) lim - / P{Qi(s) = m - 1, Q 2 (s) =j + l}ds 

(8.14) '7° 7 ° 

1 /"* 

= lim lim- / P{Q(s) = m + j + l}ds, 

and J 2 — J\ = o(l) as — > 00. The theorem is proved. □ 

Note, that the analogue of Theorem 8.2 for the Markovian multiserver 
retrial queueing system is proved in Falin, and Templeton [17]. 

9. An algorithm for numerical calculation of the model 

The aim of this section is to provide the method for calculating 

(9.1) lim \ f t F{Q 1 (s)=i,Q 2 (s)=j}ds. 
t-*oo t J 

It is worth first noting, that approximation of ()9.1|) with the aid of simulation 
by straightforward manner is not an elementary problem. Taking T large 
and approximating (|9.1f) by 

(9.2) 1 j\{Q 1 (t)=i,Q 2 (t)=j}dt 

is not realistic. For satisfactory approximation of the integral by sum it is 
necessary to take a small step A. Then number of terms should be very 
large, and the computational procedure becomes complicated. 

A more simple way is to use the Lebesgue theorem on dominated conver- 
gence: 



(9.3) 



1 [* 

lim - / P{Qi(s) =i,Q 2 (s) = j}ds 

t^oo t J Q 

1 i l 

= P lim - / I{Qi(s) = i, Q 2 (s) = j}ds. 



In this case, taking T large enough we estimate 

(9.4) i [ T I{Qi(t)=t,Q 2 (t)=j}dt 
1 Jo 

rather than (|9.2|1 . The trajectories of I{Qi(t) = i, Q 2 (t) = j} are step- wise, 
and therefore the computation procedure of 1)9.4(1 with the aid of simulation 
is much simpler than that of ((9.2() . 

On the other hand, paying attention that relations l(5.17() - ((5.19|) contain 
the terms 

(9.5) lim il f t I{Q 1 (s-)=i,Q 2 (s-)=j}dA(s), 
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then approximation of 1)9.5)1 by 



(9.6) 



^E / I{Qi(t-) = i, Q a (t-) = j}dA(t) 



is in turn simpler than approximation of (|9.4|) . 

Indeed, 1)9.5)) and (|9.6() are the Stieltjes-type integrals. That is for a given 
realization of A(t, u>) they can be represented as finite sum in the points 
of jump of A(t,uj). Then, the symbol E in (|9.6|) requires averaging of the 
results after a large number of realizations. By the Lebesgue theorem on 
dominated convergence we have 



rather than (|9,6|) . This means that it is sufficient only one long-run simu- 
lating. 

The main difference between the computation procedures for (|9.4j) and 
1)9.8)1 is the following. Whereas 1)9.8)1 requires to compute only the number 
of jumps in interval (0,T), Q9.4JI takes also into account the lengths of the 
time intervals that the process spends in states (i, j). This has no essen- 
tial significance for the algorithmic complexity of the simulation program. 
However, the time intervals that the process spends in phase states may 
vary in wide bounds, that is the average time that the process spends in 
different states (i, j) and (k, I) may have a large difference. As a result, 
1)9.4)1 is sensitive to these variations in the sense, that a small error for time 
average in specific state can result an essential error of (|9.4j) . Particularly, 
1)9.4(1 is sensitive to the behavior of the process in the boundary states (0, i), 
associated with the case where the main queue is empty. 

Hence, from the computational point of view a necessary accuracy for 
(|9.8|) can be achieved easier than that for 1)9.4)1 . Thus, between two suggested 
approaches for approximation of 1)9.1)1 . the approach based on simulating 
()9.8)l with subsequential numerical solution of the system of equations is 
preferable than a more straightforward approach based merely on simulation 



For this reason the computational procedure below is based on 1)9.8)) rather 
than on lfP|) . 

As values aij(T) are calculated, 1)5.17)1 - 1)5.19)) can be approximated as 
follows. 



(9.7) 




Therefore, for enough large T it can be taken 



(9.8) 




of 
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(i) In the case i = 

i r l 

IH lim - / P{Qi(s) = 1,0 2 W =j}ds 

l9 ' 9j - M2i Hm ~ / P{Qi(s) = 0,g 2 (s) = j}ds 

~ a j(T). 

(ii) In the case i = 1, 2, . . . , m — 1 

1 /■* 

+ 1) lim - / P{Qi(s) = i + 1, Q 2 (s) = j}ds 

t^OO t Jq 

i r* 

- Ml t lim - / P{Qi(s) = i, Q 2 (s) = j}ds 

i— oo t Jq 
(9 W) 1 /"* 

V ' 7 -^j lim - / P{0i(s) = i,Q 2 (s)=i}d S 

t^OO i Jq 

1 /•* 

+ /i 2 (i + 1) lim - / P{Qi(s = i - 1, Q 2 (s) = j + l}ds 

t-*00 t Jq 
~ a i,j( T ) ~ a i-lj(T)- 

(iii) In the case i = m 

1 /"* 

M 2 (j + 1) lim - / P{Qi(s) = m - 1, Q 2 (s) =j + l}ds 

t^OO t Jq 

{ ' )A ' ' - mm lim - f P{Qi(s) = m,Q 2 (s) = i}ds 

~ a m ,j(T) — a m _ij(T) — a mj -_i(T), 
where a mj _i(T) = 0, and a~ij(T) = 0. 

Equations l|9.9|) - l|9.11|) are similar to those for the Markovian system. The 
only difference that in the case of Markovian system the values lim^oo t^ 1 J 
P{Qi(s) = i, Q 2 (s) = j}ds are replaced by limiting stationary probabilities 
Pij (see reference (|7.4(l ). Then, the traditional way to estimate (|9.1j) is 
based on one of the known truncation methods. For example, two different 
methods are described in Chapter 2 of the book of Falin, and Templeton 
[17]. Other method can be found in Artalejo, and Pozo [6]. However, the 
analysis of the above non-Markovian system by truncation method is much 
more difficult than in the case of Markovian system. Reduction to truncated 
model implies that the initial model should be replaced by state-dependent 
model. In the case of Markovian system, the system of equations for the new 
state-dependent model is based on the Chapman-Kolmogorov equations. 
This new system of equations is an elementary generalization of the initial 
system of equations. In the case of non-Markovian model, reduction to 
truncated model leads to cumbersome analysis, and it is not clear whether 



MULTISERVER RETRIAL QUEUES 



27 



the system of equation for truncated model is similar to its variant of the 
Markovian case. 

The algorithm below provides numerical results remaining in a framework 
of the initial model. However, it is available only for the systems with 
'well-defined' parameters, when the queue-length in orbit is not large. For 
example, in the case of heavy load and low retrial rate, a queue-length in 
orbit is large, and the present method becomes unsatisfactory. 

The algorithm contains the following two steps: 
• Step 1 - initial simulation. 

The first step enables us to obtain the values a,ij(T). These values are 
then used in equations (|9.9|) - (|9.11|) . There is the finite number of equations. 
For the small e = T _1 , we define the number of equations W as W = 
max{j : aij(T) > e}. Then, according to (|9.8|) the value W is associated 
with the maximum index j for which the functional 



takes a positive integer value. For j > W the above functional is equal to 0. 
• Step 2 - solution of the equations. 

As the values ctij(T) are computed, we solve the equations and find the 
desired approximations (|9.2|) for (|9.1|) . 

Notice, that in the case where the value W is large, it is necessary to use 
one of truncation methods nevertheless. All numerical results obtained in 
the present paper are associated with the cases where W is not large, and 
we do not follow the truncation methods. 



In this section a few numerical examples for simple non-Markovian retrial 
queueing systems is provided. Specifically, the examples are provided for 
two different retrial queueing systems having two servers. One of them is 
traditional, the D/M/2 retrial queueing system. Its interarrival time is equal 
to 1. The load of the system g = (2^i) _1 varies, including low, medium 
and high load. The set of rates in the orbit varies similarly, including low, 
medium and high rates. 

The other queueing system is not traditional. Interarrival times are as- 
sumed to be correlated random variables as follows. The first interarrival 
time, £i, is uniformly distributed in interval (0,2), and the n + 1st interar- 
rival time is recurrently defined as £ n+ i = 2 — £ n , n > 1. Thus, {£ n }n>i is 
a strictly stationary and ergodic sequence of random variables having the 
uniform in (0,2) distribution, and E£ n = 1. The parameters [i\ and fj,2 of 
this system vary by the same manner as these parameters of the first system. 
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Retrial rate 


0.1 


1.0 


10.0 




Column 1 


Column 2 


Column 3 


(0,0) 


0.5988 


0.6333 


0.6343 


(0,1) 


0.0151 


0.0007 




(0,2) 


0.0017 






(1,0) 


0.3621 


0.3642 


0.3644 


(1,1) 


0.0176 


0.0001 




(2,0) 


0.0013 


0.0013 


0.0013 


(2,1) 


0.0002 







Table 1. The values of Pj j for the case of relatively low load /Ui=2.5 



The aim to consider so not standard system is the following. First, the 
system with correlated and alternatively changed interarrival times often 
appears in a large number of applications, and especially in telecommunica- 
tion networks. For example, such situation can occur when there are several 
sources, each of which sends messages by a constant deterministic interval. 
Second, our main results are related to the case of arrival point process with 
strictly stationary and ergodic increments, and it is interesting to compare 
the results obtained for this not traditional system with the corresponding 
results related to standard queue with usual, say deterministic, arrival. 

In turn, the numerical results, obtained for the D/M/2 retrial queue with 
high retrial rate, are compared with corresponding numerical results for the 
'usual' D/M/2/oo queue. The last are obtained from the known analytic 
representations (e.g. Borovkov [9]). 

For our convenience for two-server retrial queueing systems we use the 
following notation 

1 /"* 

(10.1) P id = lim - / P{Qi(s) = i, Q 2 (s) = j}ds. 

t^oo t Jq 

Note, that the limiting frequency Pjj, given by po.lj) . can be thought as 
steady-state probability. In all our numerical experiments we take T = 
100,000, and therefore e = 0.00001. 

Our numerical analysis we start from the D/M/2 retrial queueing system. 
Table 1 is related to the case of relatively low load {fix = 2.5) and different 
retrial rates. 

In the case of low retrial rate fi 2 = 0.1 by simulation we obtain W = 2. 
Specifically, ao,2(100, 000) > 0, that is the value Po,2 is positive (Po,2 ~ 
0.0017). Moreover, the maximum value of column 1 is Poo ~ 0.5988, and 
the values Poj are greater than the corresponding values Pj,-, for i = 1,2. 
This enables us to conclude that a customer, who upon arrival goes to the 
orbit, continues to spend there a long time while the main queue is empty. 
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Retrial rate 


0.1 


1.0 


10.0 


(hj) 


Column 1 


Column 2 


Column 3 


(0,0) 


0.0104 


0.1258 


0.2285 


(0,1) 


0.0155 


0.0971 


_ 


(0,2) 


0.0308 


0.0348 




(0,3) 


0.0454 


0.0155 


_ 


(0,4) 


0.0512 


0.0018 


_ 


(1,0) 


0.0104 


0.1632 


0.4675 


(1,1) 


0.0173 


0.1471 




(1,2) 


0.0322 


0.0838 




(1,3) 


0.0489 


0.0097 




(1,4) 


0.0580 


0.0041 




(2,0) 


0.0036 


0.1883 


0.2010 


(2,1) 


0.0071 


0.0642 


0.0781 


(2,2) 


0.0144 


0.0276 


0.0038 


(2,3) 


0.0211 


0.0158 


0.0008 


(2,4) 


0.0281 


0.0024 


0.0002 



Table 2. The values of Pij for the case of medium load ni=1.0 



In the case of medium retrial rate [12 = 1.0 by simulation we have only 
W = 1, that is the orbit capacity does not increases 1 at the moment of 
arrival. From column 2 of Table 1 it is seen that the values Po,i an d Pi 1 are 
sufficiently small, nevertheless Po,i > Pi,i- This can be explained by effect 
of low load. The most of time the server is empty, and the situation, when 
a customer in orbit returns to the empty queue, is typical. 

In the case of relatively high retrial rate [12 = 10.0 the simulation gives 
W = 0. In column 3 of Table 1 there are only three positive values for 
Pij which are approximately the same as steady state probabilities for the 
D/M/2/oo queueing system. 

The next table, Table 2, is associated with the case of medium load (fxi = 
1.0). In this case the value of traffic parameter g = 0.5. 

The data in column 1 of Table 2, associated with the case of relatively low 
retrial rate, show that the expected queue-length in orbit is relatively long. 
As the retrial rate increases, the expected queue- length in orbit decreases. 
Whereas for /i2 = 0.1 we have W = 18, then for fi2 = 1.0 we have W = 8 
and for /X2 = 10.0 only W = 6. 

In the next table, Table 3, the results for retrial queue with high retrial 
rate and for standard D/M/2/cc queue are compared. For the D/M/2/co 
queue it is assumed that inter arrival time is equal to 1. It is also assumed 
that \x\ = 1. Therefore the load g = 0.5. We consider the similar D/M/2/oo 
retrial queue with relatively high retrial rate \x<i = 10.0, and the comparison 
results for these two queueing systems are given in Table 3. 
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Retrial 


system 


Standard 


system 




P 


k = i + j 


Pk 


(0,0) 


0.2285 





0.2289 


(1,0) 


0.4675 


1 


0.5423 


(2,0) 


0.2010 


2 


0.1772 


(2,1) 


0.0781 


3 


0.0412 


(2,2) 


0.0038 


4 


0.0084 


(2,3) 


0.0008 


5 


0.0017 


(2,4) 


0.0002 


6 


0.0003 



Table 3. The values of Pij and Pi + j for the case of medium 
load and relatively high retrial rate 



Following the known results for the GI/M/m/oo queue given in Borovkov 
[9], Section 28, Theorem 10) and denoting 

(10.2) P = lim - / F{Q(s) = i}ds, 

t^oo t J 

we have 

(10.3) P = ^!zl, i = l,2,...,m-l, 

(10.4) Pi = ^!zi 5 i = m ,m+l, 

where A is the reciprocal of the expected interarrival time, and 

(10.5) pi = lim ¥{Q(t n -) = i}, 



n^oo 



t n is the moment of the nth jump of the point process A(t) (i.e. the moment 
of rath arrival). The explicit representation for j5j in turn can be found in 
Borovkov [9], Section 28, Theorem 9 or in Bharucha-Reid [8]. 

In our case we have: P\ = po and Pi = 0.5pj_i, i = 2,3,..., and by 
normalization condition Po=0.5(l-po). 

In turn, po = Uo — U±, pi = U\, pi = rip % ~ 2 , i = 2, 3, ip is the root of 
equation logz = 2z — 2, ip rj 0.2031, 

n i r tj r\ 1 2(1-^2) -2 
J7 = 1 - , Ui = rCi 



L \-C 2 {l-ih) 2(1-^-2 

^ = e" 1 » 0.3679, V2 = e~ 2 » 0.1353, Ci = ^i(l - V'l) -1 ~ 0.6110, 
C 2 = ^(l-^)" 1 « 0.1565, 

1 2 2(1 -^i)-l 1 2(1 - V2 



A -p Ci(l-Vi) 2(l-p)-l C 2 (l-V^) 2(1 -<p)- 2 

« 0.0823. 

Then, U » 0.8967, E/i » 0.3544. 
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Retrial 


system 


Standard 


system 




P 


k = i + j 


Pk 


(0,0) 


0.0309 





0.0455 


(1,0) 


0.1910 


1 


0.2518 


(2,0) 


0.2644 


2 


0.3101 


(2,1) 


0.1486 


3 


0.1891 


(2,2) 


0.1021 


4 


0.1334 


(2,3) 


0.0700 


5 


0.0961 


(2,4) 


0.0293 


6 


0.0311 



Table 4. The values of P{j and Pi+j for the case of relatively 
high load and relatively high retrial rate 



The first 7 values of Pk, (k = 0, 1, ...,6), are in the last column of Table 
3. In the second column of this table are corresponding values of Pi a taken 
from column 3 of Table II. The results of Table 3 are agreed with convergence 
Theorem E21 

We study now the cases of relatively high load, pL\ = 0.6. The the cases of 
high load lead to increasing queue- length in orbit, and therefore the numer- 
ical analysis becomes difficult. For example, if in addition the retrial rate is 
low, then the number of equations becomes large, and only special methods 
of analysis are necessary. For example, if \ii = 0.1, then the initial simu- 
lation shows that W > 50, and the values ao,o(100, 000), ai,o(100, 000) and 
02,0 ar e negligible. At the same time, 02,30 ~ 0.0055, 02,39 ~ 0.0331. In the 
case of medium retrial rate as \ii = 1.0 the number of equations is still large, 
W = 30. Here, the values o , (100, 000) rs 0.0014, oi, (100, 000) « 0.0049, 
o 2 ,o(100,000) ps 0.0167. The maximum value 0^(100,000) is achieved for 

1 = 2 and j = 4. Namely, 02,4 ~ 0.1265. In the case of relatively high retrial 
rate as fi2 = 10, by initial simulation we obtain W = 27. However, the 
values Ojj(100, 000) decreases in j, and the maximum value 0*^(100,000) is 
achieved for i = 2 and j = 0. Namely, O2,o(100, 000) ps 0.361. We provide 
Table 4 of some values when \i\ = 0.6 and ^2 = 10. The left side of the table 
(columns 1 and 2) contains the values for retrial queue, while the right side 
of the table (columns 3 and 4) is associated with the standard multiserver 
queue. (The results for the standard multiserver queue in this table are 
obtained by the computations analogous to that of Table 3.) 

The numerical analysis shows that, as load becomes high, the difference 
between Pjj and P%+j increases. 

Now we provide numerical results for the described above non-standard 
retrial. Recall that the first interarrival time £1 is uniformly distributed in 
(0,2). The other interarrival times are determined recurrently as £ n +i = 

2 — £ n . For this system we provide numerical example only under medium 
setting \x\ = 1 and relatively high retrial rate Hi = 10- By initial simulation 
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0,0 
0,1 
1,0 

1,1 

2,0 

2,1 
2,2 
2,3 
2,4 



Table 5. The va^ 



Not standard 
arrivals 



0.3130 
0.0880 
0.3436 
0.0420 
0.2027 
0.0100 
0.0005 
0.0001 



ues of Pij for t 



Deterministic 
arrivals 



0.2285 

0.4675 

0.2010 
0.0781 
0.0038 
0.0008 
0.0002 



le retrial systems with not 



standard and deterministic arrivals 



we obtain W = 4. This value is less than in the case of deterministic 
interarrival times (W = 6). However, whereas in the case of deterministic 
interarrival the values ao,i(100, 000) and 0x^(100,000) were negligible, in 
the case of this system a ,i(100, 000) » 0.0683 and 01,1(100,000) » 0.0382. 
In Table 5 we provide the values Pjj computed for the retrial system with 
not standard arrivals (left side of the table) and deterministic arrivals (right 
side of the table). 

The obtained results enable us to conclude, that the behavior of system 
with a not standard arrival differs from that with deterministic interarrival 
time. Specifically, in the case of the system with a not standard arrival the 
convergence of the abovementioned functionals to its limits seems slower 
than in the case of the system with deterministic interarrival time. 

11. Concluding remarks 

In this paper, an analysis of non-Markovian multiserver retrial queueing 
system is provided with the aid of the theory of martingales. The system 
of equations for this system is obtained, as well as the asymptotic analysis 
as parameter \i2 increases to infinity is provided. The representation for the 
system of equations enables us to study the system numerically, where some 
terms of the system of equation are established by simulation. Then the 
system of equation is reduced to other system of equations, similar to that 
of Markovian multiserver retrial model. The results, obtained in the pa- 
per, enable us to study not standard models of complex telecommunication 
systems arising in the real life. 
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